A double-domain spectral method for black hole excision data 
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In this paper a new double-domain spectral method to compute binary black hole excision initial 
data is presented. The method solves a system of elliptic partial differential equations in the exterior 
of two excised spheres. At the surface of these spheres, boundary conditions need to be imposed. 
As such, the method can be used to construct arbitrary initial data corresponding to binary black 
holes with specific boundary conditions at their apparent horizons. We give representative examples 
corresponding to initial data that fulfill the requirements of the quasi-stationary framework, which 
combines the thin sandwich formulation of the constraint equations with the isolated horizon con- 
ditions for black holes in quasi-equilibrium. For all examples considered, numerical solutions with 
extremely high accuracy were obtained with moderate computational effort. Moreover, the method 
proves to be applicable even when tending toward limiting cases such as large radius ratios for the 
black holes. 
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I. INTRODUCTION 

Spectral methods utilizing multiple spatial domains 
have been used by many authors in order to solve el- 
liptic partial differential equations in general relativity, 
see e.g. Q, S H S H 0, S II 13 E3 In this paper we 
demonstrate the strength of these methods by applying 
them to the calculation of initial data corresponding to 
binary black hole systems. 

In a '3+l'-splitting of space and time, Einstein's equa- 
tions may be formulated as constraint and evolution 
equations. The constraints form elliptic equations on 
spacelike hypersurfaces, which must be fulfilled by any 
data set that describes some initial state of the binary 
system. In constructing data of this kind, two differ- 
ent approaches have been explored, (i) puncture meth- 
ods [mJ3Jl3JllJF3Jlland (ii) excision techniques 

IlllTaliflElliEl^lsfllHlSl 

In the puncture methods, a special pole- like structure 
of the singularity inside the black hole is assumed which 
can be taken into account by a specific ansatz for the 
initial data. Therefore the relevant space for the con- 
straint equations is all of R 3 . In contrast, the exci- 
sion techniques solve the constraints only in the exte- 
rior of two excised spheroids within which the singu- 
larities are located. At the surface of these spheroids, 
special boundary conditions need to be imposed, see 
e.g. SHU Ed EE! Hill. A promising approach 
for providing physically realistic data that describe bi- 
nary systems in quasi-stationary orbits is given by the 
combination of the conformal thin sandwich formulation 
[2ft l27j (for a review see and the isolated horizon 
framework 000 ETEll H3- while the confor- 
mal thin sandwich equations incorporate the concept of 
quasi-stationarity into the constraint equations, the iso- 
lated horizon framework describes specific boundary con- 
ditions valid at the apparent horizons of the black holes 



that ensure a quasi-equilibrium state. 

In this paper, a new numerical scheme is presented that 
calculates binary black hole excision data by means of a 
double-domain spectral method. The scheme is an alter- 
native approach to the spectral techniques used in 0, Q 
and to those utilized in ]3 E| • It works on two spatial 
subdomains and requires relatively small spectral expan- 
sion orders to render highly accurate solutions. The com- 
plexity of the numerical scheme is comparable to the one 
presented in |16| for a single-domain puncture method, 
therefore allowing the code to run on a single processor. 

The central idea of the method is the introduction of 
two spatial domains, within which the initial data ad- 
mit a rapidly converging spectral expansion. In order 
to achieve this it is essential for the data to be smooth 
within these subdomains. At first, bispherical coordi- 
nates [13, are introduced through which the entire 
exterior of the two excised balls becomes the image of 
a single rectangular box (see Section [HJ. In particular, 
a compactification of spatial infinity is realized, which 
corresponds to a mere line on a side of the box. How- 
ever, as illustrated in Section II I II for the Laplace equa- 
tion, a solution to an elliptic equation is in general only 
C° at this line, which suggests the introduction of an 
additional mapping. This map folds the box along the 
line in question, see section llVl As a result, for a two- 
dimensional cross-section, we get a domain of pentan- 
gular shape, which we divide up into two quadrangu- 
lar ones. Each one of the two quadrangular regions is 
mapped diffeomorphically onto a square. As explained 
in Section spectral expansions for all data entries are 
carried out within these cuboids, and the collection of 
constraint equations, boundary, asymptotic fall-off, reg- 
ularity and transition conditions (the latter ones to be 
imposed at the boundary between the two domains) yield 
a complete set of equations to determine all spectral co- 
efficients. We obtain the solution by means of a Newton- 
Raphson method. For a three-dimensional code, only 
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an iterative scheme for executing the linear step inside 
the Newton solver is computationally affordable. To- 
gether with a reformulation of the regularity conditions 
(to be enforced at the axis along which the black holes 
are aligned) a specific plane relaxation scheme has been 
implemented that results in a convergent iterative proce- 
dure. In Section I VII we present first examples satisfying 
the equations and boundary conditions following from 
the above quasi-equilibrium framework. We obtain ex- 
tremely high accuracy even in limiting cases such as very 
large radius ratios for the excised spheres. Since only 
relatively little spectral resolution is needed, the code is 
of low computational cost. This makes it especially use- 
ful for the detailed study of wide classes of initial data 
sets. In Section lVlIl we discuss future applications of the 
scheme, which will include a detailed mathematical and 
physical investigation of various initial data sets and the 
use of these data for a dynamical evolution. 



II. BISPHERICAL COORDINATES 

In this section we consider bispherical coordinates 0, 
l34j through which the entire space exterior to the two 
excised spheres becomes the image of a single rectangular 
box. Let the spheres be denoted by S± and their radii 
by g± . We assume the centers of S± to be aligned along 
the x-axis at a distance D from each other and introduce 
cylindrical coordinates (x,p,ip) about the x-axis, i. e. 



y = p cos 99, z — p sin ip , 



(1) 



with (x, y, z) being cartesian coordinates. In the bispher- 
ical mapping the location of the (x = 0)-plane is chosen 
such that the product of the coordinates of inner and 
outer crossing points of the x-axis with the surfaces dS± 
(see Fig. ^| is the same for both spheres. That is 



together with 



xl-x\ = ±2g±, 



(x°, + x i ,)-(x°_ +xl) = 2D 



(2) 

(3) 
(4) 



Here clq > is the common distance of a "geometric" 
center of the two balls to the coordinate origin. It is 
given by 
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a = ^D- l ^D± - 2D*tf + + g 2 _) + (q\ - g 2 _) 2 • (5) 

The bispherical coordinates (77, £, (/j) are most easily in- 
troduced by the complex mapping 



c = iao coth 



(6) 



where 



are complex combinations of cylindrical and bispherical 
coordinates respectively. Explicitly, the coordinates p 
and x are given in terms of r\ and £ as follows: 



P 



a 



a 



sinh?] 



cosh r\ — cos £ 

sin £ 
cosh 77 — cos £ 



(8) 
(9) 



Note that the angle ip remains unchanged under this 
transformation. 

A characteristic feature of the bispherical coordinates 
is that all coordinate lines 



7/ = constant and £ = constant 



(10) 



correspond to circles in the cylindrical coordinates (x, p), 
see Fig. ^ In particular it follows that 



(x — clq coth 77) + p 2 



sinh 2 7] 



x 2 + (p- a cot£) 2 = — T7 
sin £ 



(11) 
(12) 



Among these circles we find the surfaces dS± of the ex- 
cised spheres, which are described by 



77 = 77± = iarsinn — 
Q± 



(13) 



From the preceding steps it becomes apparent that the 
entire space exterior to S± is obtained as the image of a 
single rectangular box, 



»?efo-.i7+], £e[0,7r], (^e[0,27r) 



(14) 



While we have already seen that r\ = r\± corresponds to 
the surfaces dS±, the faces £ = and £ = it (in Fig. Q] 
denoted by A± and B) represent outer and inner sections 
of the x-axis respectively. Note that spatial infinity is 
obtained as the image of the line 



77 = = £ , V e [0, 2tt) . 



(15) 



At first glance, the bispherical mapping suggests a spec- 
tral expansion of the excision data in terms of the coordi- 
nates (77, £, p>). However, as will be demonstrated in Sec- 
tion [ffi] the data are in general only C° at the above line 
corresponding to spatial infinity. We resolve this issue by 
another coordinate transformation, as will be discussed 
in Section Uvl 



III. 



BISPHERICAL SOLUTIONS OF THE 
LAPLACE-EQUATION 



We include this section about the bispherical solutions 
of the Laplace equation in order to illuminate some basic 
properties of the excision data, in particular at spatial 
infinity. Exterior to S± consider 



p + ix, C = V + i£ 



(7) 







(16) 
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FIG. 1: Illustration of the bispherical mapping 16I9H for D — 5g~, g+ = 2g- (with the azimuthal coordinate ip suppressed). 
The entire space exterior to the S± in the right panel is obtained as the image of the rectangular box displayed in the left 
panel. In both panels, solid and dashed curves correspond to constant rj- and ^-coordinate lines, respectively. The line 
corresponding to spatial infinity is emphasized by a bullet. 



subject to specific Dirichlet conditions that are enforced 
at the excision boundaries. In bispherical coordinates 
this reads as follows: 



5/2 



(cosh?y — cos£) 



CSC 2 £ + VeCOt£- -V 



where the potential V is related to U via 



U = V^/coshr] — cos£ 



(17) 
0, 

(18) 



In this formulation the solution can be found by separa- 
tion of variables 133 . It reads 



U = \J cosh r\ — cos £ x 

oo I 



;=0 m=-l 



-(l+l/2) V 



(19) 



The coefficients A; m and \xi m are obtained from an ex- 
pansion of the known Dirichlet values of V at 77 = r)± 
with respect to spherical harmonics Y^ m . 

We notice that the analytic behavior of the solution 
(11911 at spatial infinity is determined by the factor 



\J cosh i] — cos £ 



C . C 
2 sinh — sinh — 
2 2 



(20) 



V2 
2 



+e 2 +0(icn 



which is only C° at ( = 0. Although one might argue 
that the decomposition l|18(l resolves this issue for the 
Laplace equation, this is not a strategy for solving gen- 
eral nonlinear elliptic equations since the corresponding 
solutions contain terms with and without the factor (12011 . 
We therefore pursue a different approach leading to new 
coordinates in which (|20|l is analytic. 



IV. FOLDING INFINITY 

The regularity problem addressed in the previous sec- 
tion is very closely related to issues discussed in |16| . 
see formulae (48-50) therein. Indeed, introducing coor- 
dinates X and R via 



Z = R + iX = VC : 



(21) 



where the square root is taken such that X and R are 
non-negative in the computational domain (see Fig. 
immediately yields that 



\J cosh 77 — cos £ 



V2 



(X 2 + R 2 )x 



(22) 



Z 4k 



\\ j^ 4 *(2fc + l)! / 



1 + E 



Z 4k 



k=l 



A k (2k + 1)! 



is analytic in X and R, in particular at X — R — 0. 

By performing the transformation Q21JI and leaving the 
azimuthal coordinate tp unchanged, the rectangular box 
(|14|l is folded about the line l|15|) such that a pentangular 
region is created upon which X and R are defined, see 
Fig. [3 Since for the spectral method the spectral coordi- 
nates (A, B, ip) introduced below are given on a cuboid, 
we draw a border that splits the pentangular region into 
two subdomains (±) of quadrangular shape. Now the 
coordinates (A, B) with 



(A, B) e [0,1] x [0,1] 



(23) 



are mapped to values (X, R) on these subdomains via 



• region +: 

X = -X+AB + R + A sinhGu+S) + X a Be» +A 



R = -R+AB + R + Acosh(n + B) + R a Bc 



-v+A 



(24) 
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FIG. 2: Illustration of the mapping 12 U . The example dis- 
played corresponds to the geometrical picture shown in Fig. 
Q The border that separates the two quadrangular regions 
(±) is a line in this diagram. Solid and dashed curves corre- 
spond to constant A- and B-coordinate lines respectively (see 



• region 

X = 
R = 
where 



-X-AB + X-Aoasb.(jrB) + X Q Be v A 
R- AB + X-Asmh(fi- B) + R Q Bc~ u ~ A 



(25) 



Z. 



R± 



Z ± = R ± 



■ \x± = 

iX ± = ^T 1 



(26) 



fj+ = arsinh(X+/i?._) , 
fi~ = a,rsmh(R~ /X-) , 
v± = ln(X±/X ) . 



A value for Xq (Rq = tt/(2Xo)) may be chosen in order 
to create subdomains of comparable size: 



X = y/X~X+ , i? = VR^W 



(27) 



With this transformation, the surfaces dS± are obtained 
for A = 1. Exterior and interior sections of the x-axis 
(i.e. A± and B) are given by B = and B = 1 respec- 
tively, and the common border between the two subdo- 
mains is described by A = 0. 



V. THE NUMERICAL SCHEME 

Our numerical scheme possesses a number of common 
features with the method described in ^(| for the spectral 



calculation of puncture initial data. However, there are 
also specific differences, which we depict in the following 
paragraphs. 

In a spectral approximation all functions U K to be 
determined by our elliptic boundary value problem are 
considered at specific gridpoints that correspond to the 
zeros or extrema of the spectral basis functions being 
used. Here we perform Chebyshev expansions with re- 
spect to A and B and a Fourier expansion with respect 
to ip. Since the quasi-stationary framework requires so- 
phisticated boundary values at dS± for the initial data 
(see Section IVl}) . we choose to use the extrema of the 
Chebyshev polynomials so as to have gridpoints lying on 
the boundary. That is, the gridpoints are given by 



7TJ 



Bi. = sin 



2(n A - 1) 
7rfc 



2(n B - 1) 



(28) 



2irl 



where 



0<j<n A , 0<k<n Bl < I < 



(29) 



The numbers tia, ns and n v describe the spectral expan- 
sion orders of our scheme. 

The next step in the scheme is the setting of a vector 
U , whose components are derived from the values 



U jkl 



(30) 



of all functions U K at the above gridpoints (Aj, Bk,ipi) 
in the two subdomains (±). The index k 6 {1, . . . , k } 
labels the different functions appearing in the elliptic sys- 
tem, with kq being the total number of equations («o = 5 
for the quasi-stationary constraints discussed in Section 

ED- 

A way to set U that works well in the context of the 
relaxation method described below is given by the fol- 
lowing unique decomposition of U K : 



U K (A,B,<p) 
W K (A,B,Q) 



V K (A,B) 
0. 



W K (A,B,<p) 



(31) 
(32) 



We use this decomposition and collect all values V i 



K± 

and Wjfa for I > in order to build up the vector U. 
This means that in our approach the are not stored 

directly in U, but can be recovered from the entries of U 
through the above sum (JST}. 

The collection of elliptic equations valid in the exterior 
of >S±, boundary conditions imposed at A = 1 (i.e. at 
dS±), transition conditions to be imposed at A = 
(i.e. at the common border) and specific regularity condi- 
tions that need to be fulfilled at B = and B = 1 (i.e. at 
the re-axis) yield a discrete non-linear system 



f(U) = 



(33) 
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of dimension 

2K riAn B n v> . (34) 

The corresponding solution describes the desired spectral 
approximation of the solution to our elliptic boundary 
value problem. Let us discuss in detail the several entries 
fjH °f t ne vectorial function /: 

1. The spectral method enables us to calculate first 
and second derivatives of the U K from the values 
Uijt at the arj ove grid points up to the chosen ap- 
proximation order. Using them in the given elliptic 
equations evaluated at these grid points provides 
relations between the U*^. Note that the elliptic 
equations are only considered at inner grid points, 
i.e. for 

1 <j< n A -2, 

1 < k < n B -2, (35) 
< I < n v , 

thus filling 

2K {n A - 2)(n B - 2)n v (36) 
entries of /. 

2. In exactly the same manner, we use the boundary 
conditions valid at A — 1 to fill 

2n Q n B n v (37) 

entries of /. Note that the boundary conditions are 
also considered at the x-axis, i.e. at gridpoints with 
j = or j = n B - 1. 

3. Since the border between the two domains is some 
artificial transition surface that we have intro- 
duced, each function U K must be smooth there. 
In the formulation of the elliptic boundary value 
problem, this is ensured by requiring the U K to 
be continuous and to possess continuous normal 
derivatives at A = 0. These conditions arc enforced 
along the entire border except the line (|15|l where 
the values of U K at infinity are imposed, thus lead- 
ing altogether to 

2K n B n v (38) 

entries of /. 

4. Special care is needed for setting the remaining dis- 
crete equations that correspond to regularity con- 
ditions along the x-axis. Any function U that is 
smooth along the cc-axis with respect to cartesian 
coordinates can be expanded in a Taylor series: 

U(x,y,z) = U(x, 0,0) + (39) 

+p[U y (x, 0, 0) cos p + U z (x, 0, 0) sin <p] + 0{p 2 ) 



from which it follows that 

Urn U w = (40) 

and 

lim([/ p + U pw ) = . (41) 

The entire set of equations turns out to be solvable 
only if we consider both conditions 1)41) \ and l|4ip. 
In particular, we require (|40l) at gridpoints with 
<p > (i.e. I > 0) and (|41|) at gridpoints with ip = 
(i.e. 1 = 0). The corresponding relations fill the 
remaining 

4:K (n A -2)n v (42) 

entries of /. Note that different conditions need 
to be imposed if the underlying function U is not 
smooth along the x-axis with respect to cartesian 
coordinates. An example for this is given by the 
bispherical components of a vectorial function. In a 
subsequent publication we will discuss in detail this 
issue in connection with the momentum constraint 
equations. 

In treating the system l|33f) we follow very much the ap- 
proach described in |16( , see section II therein. In partic- 
ular, the solution is obtained by Newton-Raphson itera- 
tions, and the linear step inside this solver is performed 
with the preconditioned 'Biconjugate Gradient Stabilized 
(BCGSTAB)' method H^. In complete analogy to [3, 
we construct a second order finite difference representa- 
tion of the Jacobian of l|33|l . which will be used in the 
preconditioning step. Note that, because of the decom- 
position (|31|) . for every gridpoint not only adjacent neigh- 
boring points need to be considered, but also the values 
for / = 0. The resulting finite difference Jacobian matrix 
has therefore at most 27k,q non- vanishing entries per row. 

A crucial difference to the method used in [16j is the 
performance of the preconditioning step in which an ap- 
proximate solution of the system 

J FD SU = -f(U) (43) 

is obtained. The matrix Jfd appearing here repre- 
sents the finite difference approximation of the Jacobian 
J = df/dU which was mentioned above. The vector 5U 
describes a correction to U obtained through the values 
f(U)- 

We use a specific preconditioner that consists of suc- 
cessive plane relaxations with respect to the planes (p = 
constant = <pi. These plane relaxations in turn are com- 
posed of two different line relaxation schemes. We con- 
sider the method in more detail: 

1. For the first one of the two line relaxation schemes, 
let us pick one of the two regions (±) and some 
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values j, I and k. Then all values 5Wlf (or SVl* 
for I = 0) along the coordinate line 

A = Aj, B G [0,1], <p = ipi (44) 

are determined simultaneously through (|43|) while 
all other entries of 8U are held fixed. This gives a 
tridiagonal system of dimension which can be 
solved with low computational cost. 

We perform this procedure within the plane tp = 
tpi for all values j and n and in both regions (±). 
Hereby, the loop over j takes on all even values first 
and then all odd ones, thus describing a 'zebra' line 
relaxation scheme. 

2. Similarly, for the second line relaxation scheme fc, I 
and k are chosen, and all values SW^ k f (or SVj^ 
for I = 0) along the coordinate line 

Ae[0,l),B = B k ,tp = <p l (45) 

are determined simultaneously through (|43|l in 
both regions (±), thus leading to a tridiagonal sys- 
tem of dimension (2ua — 3). Here we make use of 
the continuity of the U K at the transition border. 
Note that we leave out the values at the bound- 
ary A = 1, which for specific boundary conditions 
might be problematic for the convergence of the re- 
laxation. We therefore obtain an update of these 
values only through the line relaxations described 
inHJ 

Again we apply these steps within the plane tp = tpi 
for all values k and k using a zebra loop with respect 
to k. 

3. The plane relaxation scheme within the plane tp = 
tpi is composed of several of the above orthogonal 
line relaxation schemes. A zebra loop with respect 
to I completes the relaxation method. 

The preconditioning step, consisting of a number of such 
relaxations (a typical value is 20) , is the fragile ingredient 
of the method. A simple reformulation of the boundary 
or regularity conditions might spoil the convergence of 
the scheme that otherwise is obtained through its itera- 
tive application. Since for interesting excision data the 
collection of elliptic equations and boundary conditions 
form a complicated system, it is very difficult to know 
how to pick a particular formulation in advance in or- 
der to ensure a stable convergent iteration. Nevertheless, 
for the relevant examples discussed in Section IVII it was 
possible to cast the conditions into a suitable form. 

Note that in general, relaxation schemes are compu- 
tationally quite expensive if high resolution is used. In 
this case, it would be possible to include the relaxation 
scheme within a multigrid solver. However, since the 
spectral resolution can be chosen to be rather small (see 
Section IVI|I , the method is already very efficient as it 
stands. 



A preconditiong step consisting of 20 relaxation iter- 
ations gives a good approximation for the 'BCGSTAB'- 
method. Typically, only about 6 iterations within this 
scheme are needed to complete the linear Newton step. 
Thus the numerical scheme converges rapidly to the de- 
sired finite spectral approximation of the solution. 



VI. QUASI-STATIONARY EXCISION DATA 

The numerical scheme presented is applicable to an 
arbitrary set of elliptic equations that is valid in the ex- 
terior of two spheres and subject to specific boundary 
conditions required at the surfaces of these shells. In 
particular, it should prove fruitful for the calculation of 
a variety of initial data sets corresponding to binary black 
hole systems. In this section we apply the method to the 
important example, in which the initial data are given 
through the quasi-stationary framework. 

In a first subsection we review the thin sandwich ap- 
proach which yields the set of elliptic equations to be 
considered. Next the isolated horizon boundary condi- 
tions describing a black hole in a quasi-equilibrium state 
are discussed. The main focus of this section is on the 
presentation of exemplary solutions to the corresponding 
boundary value problem, which have been obtained by 
means of the scheme to an extremely high accuracy. In 
particular, we illustrate the strength of the method in 
the limiting case of very large ratios g + / 



A. The thin sandwich equations 

In the ADM- formulation of a '3+1 '-splitting of the 
spacetime manifold, the general relativistic line element 
is written as 

ds 2 = jij (dx 1 + (Pdb) {dx j + I3 3 dt) - a 2 dt 2 , (46) 

where 7y is the 3-metric, and /3 1 and a are the shift vector 
and lapse function respectively. Einstein's field equations 
can be split up into a set of constraint and evolution 
equations for the 12 quantities 7^ and K{j where 

^ = ^(V i /9 i +V i A-ft7«) (47) 

is the extrinsic curvature (Vj represents the spatial co- 
variant derivative associated with 7y). The definition 
(|47|l yields one of two sets of six evolution equations. 

In York's 'Conformal Thin-Sandwich Decomposition' 
|2("il| the evolution of the metric between two neighboring 
slices t = constant is considered. More precisely, the evo- 
lution equations (|47|l make it possible to prescribe spe- 
cific values for the initial time derivative of the confor- 
mal metric 7™ , which is connected to 7™ via 7^ = "0 4 7ij 
where ip is a 'conformal factor'. Note that this splitting 
is unique only if we require some normalization for 7^ , 
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e.g. det(7„) 
decomposed: 



1. In addition the extrinsic curvature is 



1 



>K, K = K\ 



(48) 



o 

In order to obtain initial data corresponding to a bi- 
nary black hole in a quasi-stationary orbit, the constraint 
and evolution equations are considered in a comoving 
frame of reference in which the time derivatives of the 
metric quantities are assumed to be small initially. The 
thin sandwich decomposition allows us to set 

dt% = (49) 

to from which by virtue of (|47|l it 



in an initial slice t 
follows that 



A, 



if; 6 
2a 



with 



(L/% = V^- + VjA - -7y V fe /3 fc . 



(50) 



(51) 



In this formulation the constraint equations are given by 

(52) 



v 2 ?/> — — \tP 5 k 2 



12 



Vj(lL{3) ij - (L/3) lJ V j (In aip' 



' A-- A lJ 

4 ~ 
3 







0. 



(53) 



Note that in these formulae Vj and R are the spatial 
covariant derivative and the Ricci scalar respectively, as- 
sociated with the conformal metric 7^- . Indices are raised 
and lowered with respect to this metric. 

The Hamiltonian (|52l) and momentum constraints i|53|) 
are elliptic equations which can be used to determine the 
conformal factor ip and the shift f3 l respectively. More- 
over, in the quasi-stationary framework it is possible to 
consider an additional equation which is obtained from 
the evolution equations through the requirement 



d t K = 



(54) 



on the initial slice. This also gives an elliptic equation 
which determines the lapse function a: 

V 2 (cw/0 = (55) 



In 
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^K 2 
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' A- ■ A 1 ? 



Thus in total we have five elliptic equations (|52l 1531 155|) 

in order to determine the five quantities ip, f3 l and a. 
Note that the conformal metric 7^ and the trace K are 
free data that enter these equations. From the framework 
being presented no further restrictions on these quantities 
follow. In this paper we choose examples that correspond 
to maximal and conformally flat data, i.e. we set 

K = (56) 

and 

7<j = 5 ij (57) 

in the cartesian coordinates (x,y,z) (with 5ij denoting 
the Kronecker symbol). 



B. Boundary conditions 

The above quasi-stationary formulation is completed 
by a set of boundary conditions that control the data 
at the excision boundaries and at infinity. The excision 
boundary conditions are give n th rough th e 'Isolated Hori- 
zon Framework' [H IU1! S (see also 0, 0) 
and describe black holes in a quasi-equilibrium state. In 
particular the following is required: 

1. Within the initial slice, both excision boundaries 
are apparent horizons, i.e. two-dimensional hyper- 
surfaces with S 2 topology and the property that 
the outgoing null vectors k possess vanishing ex- 
pansion. 

2. Initially, the apparent horizon is tracked along k 
and its coordinate location does not move in the 
time evolution of the data. 

Cook 0, |2{| incorporated these requirements into the 
thin sandwich formulation and arrived at the following 
boundary conditions for ip and (3 l that are required at 
dS ± : 



S*Vi0n^) 



1 -2 ~i 

atp s ■ 



■ft 



(58) 
(59) 



Here the vector P is the outward pointing unit vector 
normal to dS± (with respect to jy) and hij — Jij — SiSj 
the conformal metric induced on dS± . The quantity J is 
given by 



ij%jA^ + -K 



(60) 



and the vector 0f, tangent to dS± is proportional to a 

conformal Killing vector of the conformal metric h lJ , with 
the proportional factor describing the angular velocity of 
the black hole. Note that the isolated horizon frame- 
work does not yield a boundary condition for the lapse 
a, which we are therefore free to choose. 

In the examples discussed below we consider corota- 
tional black holes for which the angular velocity of the 
black holes, as seen within the comoving frame of refer- 
ence, vanishes, i.e. 



Pn = . 



(61) 



Moreover, we take the following boundary condition for 
the lapse, which was used among others in [Toj : 



FViiaip) = 0. 



(62) 



At infinity we impose the appropriate asymptotic be- 
havior 



lim tj) = lim a = 1 

r I — * oc r I — > 00 

Jim [(3 1 - (fix r ) 1 } = 

\r\ — *-oo 



(63) 
(64) 
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resulting from the fact that (3 l is a shift vector in a co- 
moving frame of reference that rotates with the angular 
velocity fi. In the cartesian coordinates r = {x,y,z) we 
take 

Q = ne z (65) 
corresponding to an orbital motion of the black holes. 



C. Examples 

A well known maximal and conformally flat solution 
to the above quasi-stationary boundary value prob lem is 
given by the Misner-Lindquist initial data :23, 36] which 
are used for the evolution of a head-on collision of two 
black holes. These data are characterized by vanishing 
a and j3 1 ^ at the horizon (i.e. these data do not fulfill the 
condition i|62|) which we impose below) and moreover by 
£1 = (from which it naturally follows that they do not 
represent two orbiting black holes in a quasi-stationary 
state). In order to find quasi-stationary data, we need ad- 
ditional requirements such as the equality of ADM and 
Komar masses discussed below. Nevertheless, we start 
the investigation of the quasi-stationary boundary value 
problem in question by calculating solutions that rep- 
resent modifications to the Misner-Lindquist data. In 
particular, we take the boundary condition (|62|l (instead 
of a. = 0) but retain all remaining conditions. Whereas 
the corresponding boundary value problem for a single 
black hole admits a solution which can be given explic- 
itly |lCt l3?| . no exact solution is known in the binary 
case. However, we may take a superposition of the solu- 
tions known for the single black hole in order to create 
'start up' data for the Newton-Raphson scheme inside the 
numerical scheme. We consider the resulting solutions 
for D = 10g_ and four different choices of the radius 
ratio, q+/q- £ {1; 10; 100; 1000}. The excellent conver- 
gence of the spectral method has been checked globally 
for all functions U K appearing in our elliptic boundary 
value problem. We illustrate it in Fig. [3| for the total 
ADM-mass of the system which is given in Cartesian co- 
ordinates by the following surface integral evaluated at 
infinity: 



5(Madm) 



M a 



DM 



1 

16tt 



(7y,j - ljj,i) d2S i 



(66) 



We find that the convergence seems to be geometric, as 
exhibited by a roughly linear decrease of the error in this 
diagram. 

Note that almost machine accuracy is reached for all 
ratios considered, thus proving that the method is well 
suited to the case of extreme radius ratios (and likewise 
to situations in which the distance D is extremely large). 
This may be clarified by the following considerations: 

We know from the study of single black holes that the 
Dirichlet boundary values of the data remain restricted 
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FIG. 3: The convergence of the ADM-mass correspond- 
ing to binary black hole initial data sets with vanishing or- 
bital angular velocity Q — 0. The geometrical parameters 
of the individual configurations are given by D — 10q- and 
Q+/Q- G {1; 10; 100; 1000}. For these axisymmetric calcu- 



lations riA = ns = n and n v 



3 have been chosen. We 



compared the corresponding results for the ADM-masses to 
those of reference solutions with n G {50; 60; 70; 80} for the 
above choices of q+/q-. 



independent of the ratio Q+/ Q-. However, because of 
the different scalings involved in the problem, one might 
not expect that the normal derivatives of the data at the 
surface of the small sphere remain bounded as g+ / Q- in- 
creases. This in turn would mean that a spectral expan- 
sion needs more and more resolution, which for extreme 
ratios becomes computationally unachievable. 

To illustrate why nevertheless our method handles 
these critical cases with moderate computational effort, 
consider a special solution ifH?)l to the Laplace equation 
(|T7|l . given by 



V 



(67) 



It corresponds to the potential of a single point mass lo- 
cated at (x,p) — (— a ,0) and possesses a finite value 
V = 1 at <S_. For the potential V we need not fold the 
rectangular box l|14|) at the line (|15fl in order to obtain 
analyticity there, but may directly introduce spectral co- 
ordinates on (|14fl . i.e. (A, B) via 



V = V- + (v+ ~ f]-)A. 



£ = 7rB. 



Now the derivative 
dV 



(68) 



(69) 



A=0 



tends only logarithmically to infinity as g-/ao ~~ * (see 
(|T3jl h Since the folding of lfT^|) along (|T3)l does not mod- 
ify the qualitative behavior of the derivatives, a rapid 
spectral convergence emerges generally in these critical 
cases. Note that our coordinates are somewhat similar 
to those introduced in |3Sl ] via a logarithmic mapping in 
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order to capture very different length scales appearing 
there. 

In the second example we calculate initial data cor- 
responding to two corotational black holes in a quasi- 
stationary orbit. It has been argued |'l Ii ll(ll25l] that a 
suitable value for the angular velocity SI is obtained by 
requiring the equality of the ADM-mass and the Komar 
mass which is defined by the following surface integral at 
infinity: 

M k = ^ / Y J (Via - l3 k K lk )d 2 S . (70) 

oo 

In the example considered here, D — 10g_ and g + = 
have been chosen. For these data an orbital angular ve- 
locity SI w 0.036975/g- emerges. Note that we introduce 
the vector 

/? 4 = I3 l - (Si x r) 1 (71) 

into our numerical scheme which enables us to set definite 
values of the corresponding U K at infinity. 



<5(M A dm) 




10 15 20 25 30 35 40 



n 

FIG. 4: The convergence of the ADM-mass corresponding 
to a corotating binary black hole initial data set in a quasi- 
stationary orbit. The geometrical parameters of the configu- 
ration are given by D — 10g~ and g+ — g-. The data are 
characterized by the equality Madm = Mk through which 
an orbital angular velocity of Q w 0.036975/f?_ emerges. 
For these calculations the spectral resolutions ua = tib ~ 
2n v + 1 — n have been chosen. We compared the correspond- 
ing results for the ADM-mass to that of the reference solution 
with n — 51. 

Again we obtain a rapid overall convergence of the 
method which we illustrate by displaying the relative 
error of the total ADM-mass, see Fig. (J1J. In contrast 
to the above axisymmetric examples, a roughly linear 
section of the curve for small resolution is followed by a 



much more flattened part as the accuracy approaches 10 
digits. A similar convergence rate was found for puncture 
initial data possessing individual linear momenta of the 
black holes [lfj . We plan to discuss mathematical issues 
of this kind thoroughly in a subsequent publication. 



VII. CONCLUSIONS 

In this paper we presented a numerical scheme that 
calculates binary black hole excision data by means of 
spectral methods. The central idea of the scheme is the 
introduction of specific coordinates that are related to 
bispherical coordinates, in order to permit rapid conver- 
gence of the spectral expansions. In particular, the entire 
space exterior to the two black holes is obtained as the 
image of two spatial domains within which the spectral 
expansions are carried out. A special formulation of the 
boundary and regularity conditions enables us to use a 
specific iterative technique to approach the solution that 
corresponds to the spectral approximation of the desired 
data. 

The scheme has been used to calculate examples cor- 
responding to the important quasi-stationary framework 
which is given by the conformal thin-sandwich decompo- 
sition together with the isolated horizon boundary con- 
ditions. It exhibits a rapid spectral convergence of the 
numerical solutions up to an extremely high accuracy. In 
particular, configurations with large radius ratios of the 
black holes may be considered up to this precision, with 
only moderate computational effort. 

In future applications of the scheme we shall calculate 
a variety of initial data sets. The high accuracy of these 
data will enable us to investigate physical and mathemat- 
ical properties thoroughly. Moreover, we intend to study 
several formulations of the constraint equations in order 
to select physically interesting data which correspond to 
a binary system with a realistic content of ingoing and 
outgoing radiation. For this it will prove fruitful to han- 
dle explicitly extreme configurations that correspond to 
limiting cases (such as the test mass limit). Finally we 
plan to use the data in a dynamical evolution of the sys- 
tem which ultimately will help clear up the physical sig- 
nificance of the data being considered. 
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